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Difference imaging is a technique for obtaining precise relative photometry of variable sources in crowded stellar fields 
and, as such, constitutes a crucial part of the data reduction pipeline in surveys for microlensing events or transiting extra- 
solar planets. The Optimal Image Subtraction (OIS) algorithm of Alard & Lupton (1998) permits the accurate differencing 
of images by determining convolution kernels which, when applied to reference images with particularly good seeing and 
signal-to-noise (S/N), provide excellent matches to the point-spread functions (PSF) in other images of the time series to 
be analysed. The convolution kernels are built as linear combinations of a set of basis functions, conventionally bivariate 
Gaussians modulated by polynomials. The kernel parameters, mainly the widths and maximal degrees of the basis function 
model, must be supplied by the user. Ideally, the parameters should be matched to the PSF, pixel-sampling, and S/N of 
the data set or individual images to be analysed. We have studied the dependence of the reduction outcome as a function 
of the kernel parameters using our new implementation of OIS within the IDL-based TRIPP package. From the analysis 
of noise-free PSF simulations of both single objects and crowded fields, as well as the test images in the Isis OIS software 
package, we derive qualitative and quantitative relations between the kernel parameters and the success of the subtraction 
as a function of the PSF widths and sampling in reference and data images and compare the results to those of other 
implementations found in the literature. On the basis of these simulations, we provide recommended parameters for data 
sets with different S/N and sampling. 
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1 Introduction 

Many astrophysical experiments rely on the precise deter- 
mination of lightcurves from sources which are either weak, 
weakly variable, and/or situated in densely populated back- 
grounds. Prominent examples are the detection of extrasolar 
planets via the transit technique, gravitational microlensing, 
and supernova searches. As these events are intrinsically 
rare, dedicated large-scale surveys have only become fea- 
sible due to automated reduction and photometry of large 
amounts of CCD data. An overview of microlensing sur- 
veys (Ogle, Eros, Macho, Moa, and Planets) is given 
in Dominik et al. (2002) and the references cited therein. 

Whenever stellar point spread functions (PSFs) overlap, 
the determination of the background correction for aperture 
photometry is difficult at best and heavy blending of stel- 
lar profiles renders the assignment of flux to one source or 
another ambiguous or even impossible. Iterative deblending 
methods of different kinds have been implemented and are 
in wide use: photometry packages like Daophot (Stetson 
1987) or SExtractor (Benin & Arnouts 1996) succeed in 
obtaining a high level of precision. Nevertheless, their accu- 
racy decreases as the degree of blending increases and in the 
densest Galactic fields, it is impossible to deblend without 
additional information. 
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If the source is variable and the background is (roughly) 
constant, then changes in brightness - if not the absolute 
brightness - can be measured if one can successfully sub- 
tract the non- variable background by taking the difference 
between the images, coiTected for any differences in scale 
and seeing; such an analysis is called Difference Imaging. 
If the images of a time series are compared against a refer- 
ence image taken from the same series, the difference im- 
ages obtained by adequately subtracting the reference im- 
age should be empty of any signal except for a few variable 
objects protruding from the background as positive or neg- 
ative brightness variations. Their flux relative to the value 
defined by the reference can then be measured by more clas- 
sical aperture photometry. For difference imaging to work, 
the PSFs in the reference image and in data image thus have 
to be matched exactly. Gould (1996) and Tomaney & Crotts 
(1996) first attempted difference imaging by adapting the 
data image to a reference given by the broadest PSF. As this 
deconvolution method deteriorates image quality, it worked 
only with the highest S/N data. The nonlinear PSF fitting 
introduced by Kochanski, Tyson & Fischer (1996) was more 
robust but numerically time-consuming. 

To date, the most successful difference imaging algo- 
rithm is Optimal Image Subtraction (OIS) first sug- 
gested by Alard & Lupton (1998; hereafter AL98) and im- 
plemented in their Isis pipeline. Starting with Alard (1999), 
it has successfully been applied in many different kinds of 
photometric surveys: microlensing campaigns (Alcock et 
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al. 1999); surveys for variable stars (Olech et al. 1999); su- 
pernova searches (Mattila & Meikle 2001); transit planet 
searches (Mallen-Onalas et al. 2003); etc. 

Optimal Image Subtraction determines convolution ker- 
nels which transform reference images into the data images 
via a linear least-squares fit to a pre-defined set of basis 
functions, cleverly avoiding the problems of non-linear fits 
to the parameters of a particular PSF form. The reference 
image is either that with the best seeing (and S/N) or a 
coaddition of several good images. Because of its linearity, 
OIS permits the processing of whole images and hence uses 
all available information. 

Although the convolution kernels are linear combination 
of pre-defined basis functions, there are still "external" pa- 
rameters which have to be supplied before the fitting proce- 
dure can be started; these parameters are explained in detail 
in Section |2] We investigate the dependence of difference 
image quality on the values of these parameters in Sections 
[3] and ID using different sets of simulated data and conclude 
with an outlook in section |5] 

2 Optimal Image Subtraction 
2.1 The OIS algorithm 

Strictly speaking, the algorithm presented by AL98 does not 
match PSFs but whole images: the convolution of the refer- 
ence image R{x,y) with a suitable kernel K{u,v) results 
in a model image J{x, y) which represents the best approx- 
imation in the sense of -fitting to the data image I{x,y). 
The background differences S{x, y) between the images are 
also fitted simultaneously: 

/(x, y) = J[x, y) = y) (g> K{u, v) + S{x, y). (1) 

While OIS thus does not require any isolated stellar PSF to 
be retrievable out of the blended profiles in the images, it 
is still more accurate the more distinguishable sources the 
images contairQ. 

The convolution kernel K{x^ y) is a linear combination 
of basis functions B*^*' (u, v) with u and v denoting the PSF 
kernel coordinates: 

Wp3t-1 



=0 



(2) 



The model image is then: 



Jix,y)= J2 a^'\Rix,y)-B^'Hu,v)')+Six,y) . (3) 

i=Q 

The convolution being linear, the model images can be ex- 
pressed by the same fit parameters (a^*'') as a linear combi- 
nation of iVpsf kernel and iVbg background basis images 

J{x,y)^ J2 a(''>C^'H^,y) (4) 

i=0 



' In real life, relatively isolated stars remain indispensable for image 
registration in the existing photometry pipelines. 



with a total number of = iVpsf + iVbg free parameters. 

Inserting eq. |4] into the definition of the estimator 
and approximating the pixel count errors a{x, y) to be nor- 
mally distributed, one obtains a linear system of equations, 
the normal equations: 



c = M • a 



(5) 



with a being the vector of N fit parameters. The elements 
of the vector c and matrix M are defined as follows: 



(cT{x,y)-^ ■C^'\x,y)-C^'\x,y)) , (6) 
Y,(a{x,y)-^ ■ I{x,y) ■ C^^{x,y)) . (7) 



E 

x,y 



2.2 Basis functions according to Alard & Lupton 

There is no mandatory choice for any particular OIS ba- 
sis functions. AL98 define their basis functions B^^^x, y) 
as Gaussians of G different fixed widths {hg) multiplied by 
polynomials in kernel coordinates: 



B^\u,v) = exp I - 



2hl 



(8) 



For the exponents j and fc, the relations < j,k < dg 
and < j + fc < dg hold, where dg is the maximal de- 
gree of the g-th basic Gaussian component. The multi-index 
i = {g,j,k) comprises the indices running over the basic 
Gaussians and the polynomial exponents. The widths bg de- 
termine how much a PSF's width increases upon convo- 
lution with the respective kernel and are therefore called 
broadening parameters. Together with the maximal degrees 
dg, the bg are the external parameters to be supplied before- 
hand. Further adjustments can be made to the number G of 
Gaussian components and the size Al of the array represent- 
ing the convolution kernel for computation. The background 
difference between reference and data image is expanded 
into a polynomial in pixel coordinates with a maximal de- 
gree dbg defined the same way as the (dg). Thus, the basis 
images are the following: 



psf 
psf 



(9) 



In the framework of an improved algorithm, allowing 
for a spatial variation of the convolution kernel over the chip 
to model PSF variations, Alard (2000) redefined the basis 
functions by subtracting the z = function from every other 
basis function of nonzero integral. Our difference imaging 
pipeline within the IDL Tripp package (see below) follows 
the latter definition of the basis functions. Further indepen- 
dent implementations of OIS in the literature share the def- 
inition of basis functions from AL98 or Alard (2000): e.g. 
Wozniak (2000; hereafter WOO); Bond et al. (2001); Gossl 
& Riffeser (2002; hereafter GR02). 
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2.3 Difference Imaging in the Tripp package 

The Time Resolved Image Photometry Package, Tripp, is 
an IDL data reduction package for the automated process- 
ing of large CCD time series. Up to now, it has mainly been 
used for differential photometry with a clearly defined ob- 
ject of interest. Difference image analysis, by focusing on 
the detection of variable sources in crowded (and mostly 
larger) fields, brings somewhat complementary specifica- 
tions into play. We have added to the Tripp pipeline, as 
described in Schuh et al. (2003), an alternative branch of 
data flow. New top-level routines have been added for the 
interpolation of all images to a common coordinate grid, 
the actual PSF matching and image subtraction as well as 
image coaddition. Shared information, e.g. the external pa- 
rameters, is provided by an adaptation of the log file control 
from the original pipeline. 

Before resampling an image to the reference grid, sat- 
urated pixels and pixels too close to an image edge to use 
them for convolution are detected and stored in a bad pixel 
mask. Pixels in the new grid to which those bad pixels in 
the original grid contribute are flagged in the final bad pixel 
mask and left out in the computation of the function. 

The sums in eqs. (|6]l and (|7]i are evaluated on rectangu- 
lar subframes of the data images, similar to the methods of 
AL98 and Alard (2000). In the normal (fast) mode, one sin- 
gle convolution kernel is obtained from the combined data 
in all subframes. In order to account for spatial PSF varia- 
tion or regions of particular interest, the user can choose to 
determine "local" kernels for some or all of the subframes 
which then rely only on local PSF information. 

Variable sources are identified by running an adapta- 
tion of the Daophot find function (Stetson 1987) on a 
weighted sum of difference images. For these significantly 
variable sources, flux differences to the reference frame are 
measured by Tripp aperture photometry of the difference 
images. The flux scale in the difference images is that of the 
reference frame; thus, difference lightcurves may be cali- 
brated by aperture photometry of the reference. 

The following investigations on the outcome of differ- 
ence imaging runs using different sets of external parame- 
ters have been carried out using the Tripp difference imag- 
ing pipeline. In some analyses, detailed below, simulated 
images from the Alard & Lupton Isis package were pro- 
cessed. The consistency between results from Tripp and the 
ISIS pipeline (written in C) was assessed by running both 
codes with the same parameters over those images. Except 
for the more elaborate handling of edges in Tripp, the re- 
sulting difference images are similar to the level of showing 
the same noise pattern from residuals of source subtraction. 

3 The vectors of maximal degrees 

3.1 The maximal degrees and number of parameters 

The number of free parameters is determined by the maxi- 
mal degrees of the modifying polynomials for the G Gaus- 



sian components of the convolution kernel. With A(n) ~ 
Sfc=i being the n-th triangular number, the number of pa- 
rameters associated with a vector (di, . . . , c?g) of maximal 
degrees is ^{dg + !)• (This can be verified by count- 

ing the possible combinations for j and k in eq. ([8]l.) The 
simultaneous fit of the differential background adds further 
parameters. The benefit in x^^^ of using higher numbers N 
of basis functions must be traded against the runtime of the 
difference imaging code. The N x N normal equation ma- 
trix and iV-element vector representing the linear problem 
of x'^-fitting lead to a quadratic increase of computations 
with parameter number Thus, the relation between the Xicd 
and n is of practical importance. 

3.2 The simulated point spread functions 

The crucial element of OIS is matching one PSF onto the 
other by applying an appropriate convolution kernel. In ab- 
sence of intrinsic variability, a data image containing an 
isolated stellar profile should ideally be reproduced by the 
model image resulting from PSF matching such that the dif- 
ference image will show zero flux. The fit's therefore is 
a direct measure of the residuals from PSF matching. 

Many of our tests on the OIS algorithm have been car- 
ried out using simulated test data of single objects. Both 
reference and data images for each set contain a single PSF 
without noise in a 64 x 64 pixel frame. Because our inves- 
tigations are of the influence of external parameters and not 
the quality of OIS on realistic data per se, neither intensity 
and position offsets between images nor noise have been 
introduced in those data. 

While the reference PSF was chosen to be pure Gaus- 
sian by construction, the data profile function consists of a 
weighted sum of a Gaussian and a Lorentzian function with 
the same center position: in terms of its radial distance r, it 
is given by 

<^{r) = (l-O-exp (-^) + ^ (10) 
V ^SqJ l + (r/sL) 

with sq and sl being the width scales of the Gaussian and 
Lorentzian component, respectively. Throughout this arti- 
cle, the width of a profile function will be defined as the geo- 
metric mean of the standard deviations s of fitting Gaussian 
exp (r^/2s^) along the Cartesian axes. Adopting ( — 0.2 is 
a reasonable representation of a seeing profile with a cen- 
tral core and weak extended wings; we will refer to this 
PSF as LORENTZ20. The LORENTZ20 widths were cho- 
sen to be SG — l-2srcf and sl = l-5srcf, with Siof the 
width of the Gaussian reference peak. This yields a relation 
Sdat = l-22S].ef between the final widths in both images. 

Due to its Lorentzian shape at large centroid distances r, 
an isolated LORENTZ20 PSF can be measured at radii where 
it would be completely dominated by noise under realis- 
tic circumstances. To localise the PSF, we multiplied the 
LORENTZ20 PSF with an Gaussian envelope function in 
noise-free simulations. It confines the PSF to the size of the 
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image displaceme 



Fig. 1 Dependence of the reduced on displacements 
between frames, i.e. residuals of improper image registra- 
tion. The curves represent the results for different sam- 
plings using the LORENTZ20 crowded field simulations dis- 
cussed in sec. 14.31 The sampling in the data frame was 
Sdat = 1-95 px in the continuous curve, Sdat = 2.85 px in the 
dashed curve, and s^at = 3.6 px in the dash-dotted one while 
the sampling is Sicf = 1.8 px in the reference frame. 



4 - 



AL: d = (6,4,^) 
W: d = ('l,3,2) 
GR: d = (6,4,2) 



AL 



a: d = (4,2,2) 

b: d = (4,4,0) 

c: d = (2,2.2,0) 

d: d = (2,2,2,2) 

c: d = (4,2,2,0) 



20 30 40 50 
number of free parameters 



70 



Fig. 2 The minimal reduced for the LORENTZ20 data 
in the parameter range 0.5<6i<4.0 pixels as a function of 
the number N of free parameters given by the maximal de- 
grees assigned to the Gaussian kernel components of differ- 
ent width: G = 3 (crosses); and = 4 (diamonds). Results 
for the triplets d = (4, 3, 2) (WOO), d = (6, 4, 2) (GR02) 
and d — (6, 4, 3) (AL98) are represented by the initials of 
the respective authors. 



kernel without changing the relevant central part of the pro- 
file. 



3.3 Tlie number of free parameters 

In reductions of perfectly aligned test data, basis function 
with an odd number either for j or k contribute nothing to 
the kernel solution. In real data, these odd basis functions 
permit to process irregular or asymmetric PSFs or can shift 
the PSF centroid and thus compensate for any residual mis- 



alignment between images. Test reductions for which an ar- 
tificial displacement between reference and data frame had 
been introduced showed little increase in x^cd misaUgn- 
ments < 1 px compared to the case of perfect image regis- 
tration (fig. [T]i3 Increasing a maximal degree from even to 
odd is found to decrease Xrcd ^^^^ '^h™ increasing in from 
odd to even. For that reason, we will only consider even dg. 

In the implementations discussed in literature (AL98; 
WOO; GR02), the number of Gaussian components is G = 3; 
featuring decreasing maximal degrees dg with increasing 
broadening parameters bg . Kernel components having high 
polynomial exponents are larger in extent than basis func- 
tions with lower exponents while the most important part of 
the PSF model is the center. Therefore, in our study, the dg 
decrease or remain constant with increasing bg. 

We tested the effects of diverse vectors of maximal de- 
grees (dg) using Tripp. Figure |2] shows the minima of x^ed 
found in the interval 0.5 px < 6i < 4.0 px as a function of 
the number N of free parameters for all vectors of G = 3 or 
G = 4. Results from similar runs using the vectors suggested 
by AL98, WOO, and GR02 have been added and marked by 
the authors' initials H 

The most noticeable feature in fig. |2] are three points at 
significantly higher x^cd '■h™ most of the others. These be- 
long to d= (2, 2, 0) at A^=13, d=(2, 2, 2) at A^=18, and 
d=(4, 2, 0)atA^ = 22. Note that the vector d= (2, 2, 2, 0) 
(case c in fig. |2]l yields a much better Xrod ^ ~ 
a rule of thumb, at least 20 parameters are needed for suffi- 
cient complexity of the set of basis function to match point 
spread functions. Beyond this threshold, x^d rather slowly 
improves with increasing n. Between N Ki30 and N 50, 
the use of four widths seems to yield better results than hav- 
ing G = 3, while there is no difference for larger N. 

The efficiency of a certain choice for (dg) can be esti- 
mated considering Xred/^^- maximal values for this 
quantity are obtained for those vectors denoted a to e in 
fig. ID This selection also shows that extending (dg) by a 
maximal degree of zero (a pure Gaussian basis function) 
will improve Xred' demonstrated by comparing the solu- 
tions d = (2, 2, 2, 0) vs. d = (2, 2, 0) with just one addi- 
tional parameter 

The five most efficient vectors have maximal degrees 
dg < 4. This corresponds to j = 2, k — 2 being the last qual- 
itatively different distribution of flux enhancements and de- 
pressions in the convolution kernel to be added with increas- 
ing maximal degree. Tliis implies that a further reduction 
in the PSF matching residuals via a higher number of 
fitted parameters can be most easily obtained by using 
more Gaussian widths and maximal degrees up to four. 
In the rest of this article, we will investigate d— (4, 2, 2) in 
greater detail. 



This means that one could - in principle - skip the image registration 
if the offset between images is sufficiently small and constant. 

' Due to the noise-free nature of LORENTZ20 data, the absolute scale of 
X?od defined, but can still be used as a relative measure of goodness 
as a function of sample size and parameter number 
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Fig. 3 Reduced for d = (4, 2, 2) as a function of minimum Gaussian kernel width bi obtained for LORENTZ20 data, 
kernel spacing parameter /3 = 1.5, 2.0, 2.5, 3.0, and reference PSF's widths s^ci increasing in steps of 0.1 px from s^ef = 
1.2 px (uppermost line) to 4.0 px (lowermost at 6i = 0.5 px). The corresponding data PSF widths runs from Sdat = 1.44 px 
(uppermost) to 4.28 px (lowermost), with step sizes ranging from 0.22 to 0.08 px. The convolution kernel size was fixed at 
m—19 px. Asterisks: positions of the local minima at highest bi for each sampling. Diamonds: approximate values for bi 
using 6*= 1, 6? = as in eq. (fT2] i. 



4 The kernel widths 

4.1 Kernel widths and PSF sampling 

If the number of widths of the kernel components is G = 3 
or larger, the parameter space of broadening parameters - 
not to mention the whole space for arbitrary (dg) - becomes 
complicated. The different values for the G widths (bg) used 
in the literature have been arbitrarily selected by the differ- 
ent authors, based on their own experience. 

In order to test their selection systematically (and in an- 
ticipation of the main result of this section result that Xrcd 
does not depend too sensitively on the {bg)), we adopted 
the following system: instead of choosing G broadenings 
independently, we assume they are related by a geometrical 
series such that 

bg = b,-l3a-\ 5 = 1,...,G. (11) 

There remain bi, the minimum kernel width, and /3, the ker- 
nel spacing parameter, as independent external parameters. 



The choice of a geometrical series can again be justified by 
the approximately Gaussian nature of a PSF with a steep 
intensity gradient near its centroid and wings that fade out 
super-exponentially into the background. 

If the PSF's were perfect Gaussians of widths Sicf in the 
reference image and Sdat > Sref in the data peak, the con- 
volution kernel mapping the reference onto the data would 
also be a Gaussian of width 

b - V^dat - S?ef • (12) 

This simple relation can be used as a "first-guess" value of 
bi. Given G> 1 Gaussian components and the LORENTZ20 
PSF the dependence of the optimal value for bi using the 
geometrical spacing (eq. [TTJ is more complex than the sim- 
ple relation in (eq. [T2] |. For four kernel spacing parameters 
between f3 = 1.5 and f3 = 3.0, values of Xrcd against bi are 
presented in fig. [3] the graphs show how the locations of 
local minima in Xrod niove with changes in PSF width for 
widths of the reference PSFs from Sicf = 1.2 px to 4.0 px in 
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Fig. 4 Reduced as a function of minimum Gaussian kernel width bi obtained from reduction of LORENTZ20 PSFs at 
four fixed reference image samplings s^ef — 1-5 px, 2.0 px, 2.5 px, and 3.0 px. The respective panels show Xrad(^i) f^i" ^ 
subset of the data images from fig.[3]with samplings Sdat ~ Srcf ■ A fiducial value of Sdat is given in the plots along with the 
average difference between consecutive data image sampUngs. The kernel spacing parameter is fixed at P — 2.0 while the 
maximal degrees and kernel size are set to o? = (4, 2, 2) and M = 19 px (as in fig.O. 



steps of 0. 1 px and associated data PSFs from Sdat = 1-44 px 
to 4.28 px. As a consequence of limiting the peak's extent 
by a Gaussian envelope function of constant width, there is 
no fixed relation between Sicf and Sdat and the steps in Sdat 
decrease towards wider peaks. 

Keeping in mind that PSF shapes look more Gaussian at 
high Sdat in this data set, one can evaluate general proper- 
ties of the Xrod(^i) curves. Smaller PSF widths correspond 
to seeing conditions more favourable for precise photom- 
etry. On the other hand, on a CCD with its fixed physi- 
cal pixel size, these peaks are poorly sampled by their dis- 
crete measurements. Sparsely sampled peaks impose harder 
constraints on PSF fitting than well-sampled ones. With in- 
creasing PSF width, the convex shape of the {bi ) curve 
has three minima. The global minima in x^^d' generally lo- 
cated towards the highest bi and denoted by asterisks in 
fig. [3] can lead to a drastic shift in the optimal value for 
b°^^ over a narrow range in Sdat- At higher samplings, b^^^, 
reaches its maximum and reapproaches the secondary min- 
imum located near 6i = 1.0 px nearly independent of Sdat- 
This effect is not by the changes in data PSF shape and 



width evoked by the envelope function: first guess values 
for b°^^ derived from eq. ( fT2] l and marked by diamonds only 
start to turn at much higher Sdat near the bottom of the plots. 

Comparing the results for the chosen values of /3, we 
find the most pronounced differences at low samplings. For 
increasing f3, the Xrcd 'galley at Sdat = l-44px widens to- 
wards lesser values of bi and arrives to be flat for /? = 
3.0. The observation that larger /3 correspond to lower op- 
timal values for bi and tend to result in lower gradients of 
Xicd {bi ) can easily be understood as an effect from higher 
values for 62 and 63 accompanying a given bi. This can also 
be seen from sudden increases of x^cd ^i'^h 61 which occur 
whenever a maximum of a basis function contributing rele- 
vantly to the kernel solution starts to fall onto an edge of the 
kernel array upon incrementing 61. Such function will then 
start to act as an additional fitting function to the differen- 
tial background but no longer to the actual PSF matching. 
Figure [3] shows that edges in Xrcd(^i) appear at smaller 61 
the higher the value of (3, especially for intermediate PSF 
widths. We find the choice of [3 to be uncritical with /3 > 2.0, 
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with higher (3 allowing a more convenient reduction of nar- 
row PSFs. 

Figure|4]shows the Xrcd from similar reductions at fixed 
samplings in the reference frame. The overall Xrcd optimum 
is found at a data sampling sjat slightly larger than s^ef- 
Note the ability of the basis system to build image sharpen- 
ing kernels yielding a more peaked PSF. Our results indicate 
that a 20% reduction of PSF width is possible at a reason- 
able Xred' ^ feature frequently used for PSF matching of 
good seeing frames which will be stacked to produce a less 
noisy "super-reference" for subsequent difference imaging. 

As expected, the location of the minimum moves to- 
wards higher hi with increasing Sdat at a given s^cf- In ad- 
dition, we find optimal minimum widths = 0.8 px at 
Sicf = 1.5 px increasing to 6^"^* — 2.2 px at s-^d = 3.0 px. 
Shape and width of the x^^^ valley both depend on sampling 
of the reference frame and the sampling difference between 
reference and data image. 

The detailed structure of the course of x^{hi), espe- 
cially h°^'^ , depends on the individual images to be matched. 
Therefore, after a comparison of several setups, we only dis- 
cuss the more robust qualitative features here. It should be 
mentioned, that changing the kernel size M not only re- 
duces but alters the curves for a given Sdat to a flatter 
shape with respect to bi which resembles that found for 
broader peakfl We adopted A/= 19 from AL98. 

In further tests including image registration by interpo- 
lation, the minimum value of x^od drastically increases for 
data samplings Sdat < 1 .2 px, relatively independent of Sief . 
Interpolation errors dominate the reduction outcome, thus 
limiting the utility of OIS for these undersampled data. 

4.2 Crowded field tests using isis2 data 

In order to assess the parameter dependence of difference 
imaging quality in a more realistic (and easily reproducible) 
form, we utilised a series of 50 simulated crowded field ob- 
servations contained within the ISIS package as test images. 
These data, hereafter referred to as the ISIS2 set, show a 
nearly linear decrease in seeing with index spanning a range 
from s = 1.8px to s = 3.2px. We took the first image as the 
reference and every fourth of the further ones as data im- 
ages. Averaging Xrod subseries of the latter twelve 
images, the difference imaging quality for different choices 
of bi and (3 are presented in fig.|5] 

The XTOd(^i) curves show little variation with the tested 
(3 except for bi < 0.8 px, where there is a steep growth in 
^rcd ^'^^ P ~ ^ lesser extent, for /? = 2.0 as all 

three broadenings get small relative to the PSF difference. 
With an optimal value around hi « 1.0 px, parameters in the 
interval 0.7 px < 6i < 1.6 px yield similarly good results. 
For the higher three values of /3 a secondary minimum at a 
smaller bi exists. 

To properly sample the kernel's center element an odd number should 
be chosen for M. 
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Fig. 5 Reduced as a function of hi for a simulated 
time series consisting of every fourth image in the ISIS2 
data set and maximal degrees d = (4, 2, 2). The aver- 
age over the twelve non-reference images in the series is 
denoted by (Xrcd)- Curves resulting from Tripp difference 
imaging runs with four different f3 factors as well and 
broadening parameters 6 = (0.7, 2.0, 4.0) (AL98), b = 
(0.78, 1.35, 2.34) (WOO), and b = (1.0, 3.0, 9.0) (GR02) 
are plotted. 

Overplotted in fig. |5] are results from Tripp using d = 
(4, 2, 2) with broadenings taken from the literature: b = 
(0.7, 2.0, 4.0) (AL98); 5= (0.78, 1.35, 2.34) (WOO); and 
6= (1.0, 3.0, 9.0) (GR02). By construction, the latter falls 
on the P — 3.0 curve while WOO employs a geometric spac- 
ing of /3 w V3- The AL98 parameters do not follow a mul- 
tiplicative relation, but the corresponding x^cd value is very 
close to the ones for those /3 yielding a comparable vector. 
All of these models come close to the minimum found in 
our simulations with AL98's bi — 0.7 px marking the lower 
limit for advisable bi . 

With a wider x^od valley and its minimum at b^^^ « 2.5, 
the ISIS2 reductions differ significandy from the LORENTZ 
results for a width of the data PSF of Sdat ~ 3.0 px which 
is found in the poorer sampled ISIS2 images. On the other 
hand, the agreement to the LORENTZ20 model at Sjof = 
2.0 px and f3 = 2.0 is fairly good (see fig. ^ . This can be 
explained by the fact that, in current implementations of dif- 
ference imaging, one set of parameters is chosen for reduc- 
ing the whole time series comprising images obtained un- 
der varying observational circumstances or even collected 
from several telescopes (as in collaborative microlensing 
surveys). Given that situation, the averaging selects for the 
(bg) equalising the subtraction residuals over the full ob- 
served PSF range. 

A closer inspection of the data shows that the overall 
course of the averaged Xrod t^ determined by the images 
with best seeing next to the reference image. Their reduc- 
tion requires slim convolution kernels and thus small bi. In 
order to obtain the best possible lightcurve, a low Xrcd 
especially desirable for images of high quality. This means 
applying a small bi to the whole data set. Small values of 
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61 ~ 1.0 px may not be optimal for wide data PSF relative 
to the reference, but the loss in terms of x^d compared to 
the minimum is not large. Therefore, the (bg) can be deter- 
mined by optimising them for the most sensitive data image 
of best seeing. As a positive side effect, the success of a 
certain set of (bg) can be deduced from a test run on only 
two images: the reference image and a further good seeing 
exposure. 

4.3 Crowded field tests using lorentz20 data 

Divergent dependences on 61 for ISIS2 and LORENTZ20 
data may partly be caused by the obvious differences be- 
tween LORENTZ20 and ISIS2 profiles: ISIS2 PSFs are of 
purely Gaussian shape and are highly blended, while the 
unblended PSF of our test images have extended Lorentzian 
wings, cut off by an Gaussian envelope. 

We investigated the relative importance of these effects 
in reductions of a simulated time series of images featuring 
LORENTZ20 PSFs at about the same level of crowding as in 
ISIS2. Three series of 1 + 12 images at different S/N ratios 
were created using Sdat = 1.8 px in the reference images and 
an increment of +0.15 px with each higher index in a series. 
Thus, roughly the same seeing range as in ISIS2 is covered. 
The average S/N integrated over an aperture centered on a 
relatively isolated star was set at 220 in the highest quality 
case and reduced to 100 and 30 for the other ones. 

Figure |6] shows how Xrcd depends on bi for the highest 
and lowest seeing images out of the twelve data images and 
the averages over the whole series for the S/N = 220 and 
S/N = 100 cases. Values of x^cd were normalised to the 
minimum for the specific S/N and sample subset to be able 
to present the results for different grades of seeing in the 
same plot. The S/N = 30 curves resemble their counterparts 
of intermediate quality with features smoothed out by the 
larger random component. 

In the high S/N data set, the Xicd averaged over the 
complete sample follows the best seeing data image curve 
and with 0.7 px < 61 < 1.8 px displays a similar interval 
of nearly optimal minimum kernel widths as their (3 = 2.0 
and (3 = 3.0 counterparts from the ISIS2 series (fig.©. The 
shape of the curves Xrod (^1 ) fo'" '^^e best and worst seeing 
images reproduces quite nicely the results found at both f3 
for isolated LORENTZ20 peaks of similar reference and data 
sampling. In particular, for the best seeing results at /? = 2.0, 
the positions of the minima agree quite well with the ones 
for the isolated Sdat = 2.0 px and Sicf = 2.02 px (upper right 
panel in fig.lUl. 

Comparing the heavily blended ISIS2 and LORENTZ20 
data allows us to estimate the influence of PSF shape on the 
difference imaging parameters. While pure Gaussian PSFs 
yield a relatively simple dependence Xrcd(^i)' '■^^ multiple 
minima found with LORENTZ20 data seem related to their 
more complex PSF structure. Nevertheless, both data sets 
have in common the range of usable bi . Blending itself ap- 
pears to have only a minor influence on the choice of 61. 



Table 1 Summary of the recommended settings for OIS 
external parameters. 
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With increasing level of noise, the optimum in x^^^ for 
the best seeing data image is found at larger &i, while for 
the poorest seeing image it becomes less pronounced. These 
tendencies are probably caused by the less clearly defined 
PSF in noisy images. It is unclear whether differences be- 
tween the levels for the same image at different (3 in the 
right panel of fig.|6]represent an artefact of the simulation. 

5 Conclusion and outlook 

Optimal image subtraction following the AL98 setup re- 
quires about n > 20 parameters for successful subtraction of 
constant sources (fig.|2]l. The improvement with increasing 
number of parameters is marginal; the most efficient choices 
for {dg ) are given in sec. 13. 31 A larger number G of principal 
Gaussian should be favoured over maximal degrees dg > 4. 

The choice of kernel widths mainly depends on the dif- 
ferences in PSF widths between reference and data images, 
with eq. ( fT2b giving a crude estimate. Given a sampling of 
the data image Sdat ^1-5 px, the x^d °f '■^^ reduction does 
not depend critically on the (bg). Multiplicative spacing of 
kernel widths proves to be useful, with 1 .5 < (3 < 3.0 nearly 
equally recommendable. Due to a small range in widths, 
values of /3 < 1.5 should not be used. For small differences 
between reference and data samples, high values of (3 yield 
a wider and thus more comfortable interval of low Xred- 

Difference imaging has, up to now, applied a single set 
of kernel widths to a whole time series, comprising expo- 
sures of different seeing. Parameters should then be chosen 
to achieve preferentially good reductions for high quality 
images showing a small sampling difference to the refer- 
ence. As outlined in sec. 14.21 this demands for minimum 
kernel widths in the interval 0.7 px < 61 < 1.2 px able 
to model minute PSF differences. At even smaller bi, ba- 
sis function themselves become undersampled leading to a 
higher Xred- ^ guideline for the user of OIS, these results 



© 0000 WILEY- VCH Verlag GmbH & Co. KGaA, Weinheim 



www.an-joumal.org 



Astron. Nachr. / AN (0000) 



9 




Fig. 6 Reduced as a function of S/N and minimum kernel width bi for a simulated time series of highly blended 
LORENTZ20 PSF's. Shown are Xicd averaged over twelve images with linearly increasing PSF widths (solid curves) as 
well as for the best (long dashes) and worst (dash-dot) seeing data images. Note that, with decreasing S/N, the differences 
between reductions of good and poor sampling images are also reduced. Adding further noise increasingly flattens out the 
X?cd(^i) curves. 



are summarised in table [T] One can, in principle, go on to 
define customised (bg) for individual images or groups of 
images to further optimise Xicd- 

Efforts towards this direction are likely to be impeded by 
the difficulty of deriving quantitative results on the effects 
of the parameters. The deeper reason is that the AL98 basis 
functions produce elements in the matrix M and the vector 
c from eqs. ^ and ^ typically spanning > 5 decades, in- 
dependent of the input images. Although Tripp makes use 
of the advantages of singular value decomposition, small 
changes in the image data may in some cases lead to signif- 
icantly different Xrcd values. This instability coexists with 
the general features presented above. Because this problem 
arises from large numbers introduced by the multiplication 
with polynomials of pixel indices in the ad hoc definitions 
in (eq. |9]l, it should be fruitful to look for alternative def- 
initions of basis functions. Orthogonal functions - e.g. the 
decomposition of the PSF into shapelets (Refregier 2003) - 
might have the double advantages of minimising the num- 
ber of parameters needed and simultaneously allowing for 
a wider class of PSF corrections to be applied by optimal 



image subtraction. The matching of seeing conditions with 
its subtleties continues to be a nontrivial problem which so- 
lution holds many scientific and numerical insights. 
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